Numerical modeling of quasiplanar giant water waves 
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In this work we present a further analytical development and a numerical implementation of 
the recently suggested theoretical model for highly nonlinear potential long-crested water waves, 
where weak three-dimensional effects are included as small corrections to exact two-dimensional 
equations written in the conformal variables [V.P. Ruban, Phys. Rev. E 71, 055303(R) (2005)]. 
Numerical experiments based on this theory describe the spontaneous formation of a single weakly 
three-dimensional large-amplitude wave (alternatively called freak, killer, rogue or giant wave) on 
the deep water. 
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I. INTRODUCTION 



Rogue waves are extremely high, steep and dangerous 
individual waves which sometimes appear suddenly on a 
sea surface among relatively low waves (see, for instance, 
the recent works QIHEI]. and references therein). The 
crest of a rogue wave can be three or even more times 
higher than the crests of neighboring waves. Different 
physical mechanisms contribute to the rogue wave phe- 
nomenon: dispersion enhancement, geometrical focusing, 
wave-current interaction [l|, but the most important at 
the final stage is the nonlinear self-focusing mechanism 
resulting in accumulation of the wave energy and momen- 
tum on the scale of a single wavelength Q. For weakly 
modulated periodic planar Stokes waves this mechanism 
leads to the well-known Benjamin- Feir instability 0,H|, 
generated by four-wave nonlinear resonant interactions 
2 — * 2. This instability is predominantly two-dimensional 
(2D), and it is dominant for low amplitudes (h/X < 0.09 
where h is the peak-to-trough height and A is the length 
of the Stokes wave 6]). For larger steepness parameter 
h/X, another, genuinely three-dimensional (3D), instabil- 
ity becomes dominant, which is generated by five- wave 
interactions 3^2, and results in the well-known cres- 
cent, or "horse-shoe" wave patterns (see 0, 0, E3], and 
references therein). It is important that real ocean giant 
waves are observed in situations when this 3D instabil- 
ity is not principal, and all the waves are typically long- 
crested, corresponding to a narrow-angle Fourier spec- 
trum. Thus, many essential features of a rogue wave 
formation can be observed already in purely 2D geome- 
try, as in the works by Zakharov and co-workers |2j, |3j . 
For instance, in Ref. [3j a numerical giant wave was com- 
puted with the impressive spatial resolution of up to 2-10 6 
points. Zakharov and co-workers simulated an exact sys- 
tem of dynamic equations for 2D free-surface inviscid po- 
tential flows, written in terms of the so called confor- 
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mal variables, which make the free boundary effectively 
flat (the corresponding exact 2D theory is described in 
Refs. [l(J IH El El Q EE El). With these variables, 
highly nonlinear equations of motion for planar water 
waves are represented in an exact and compact form con- 
taining integral operators diagonal in the Fourier repre- 
sentation. Such integro-differential equations are easy to 
treat numerically with modern libraries for the discrete 
fast Fourier transform (FFT) as, for example, FFTW 
|17| . Recently, by introducing an additional conformal 
mapping, the exact 2D conformal description has been 
generalized to non-uniform and time-dependent bottom 
profiles, so that a very accurate 2D modeling of near- 
coastal waves and tsunami-like processes has been possi- 
ble EE!- 

However, real sea waves are never ideally planar, and 
the second horizontal dimension might play an impor- 
tant role in the wave dynamics. Various numerical meth- 
ods have been dev elop ed for nonlinear 3D surface gravity 
waves (see 0, |2jJ, El f° r a review) . Some of them are 
based on exact formulation of the problem (the bound- 
ary integral method and its modifications; see |22l l23l |24| 
and references therein), another approach uses approxi- 
mate equations of motion, as the Boussinesq-type mod- 
els [2^, 12^1. the equations derived by Matsuno [27| . 
by Choi [28| . the weakly nonlinear Zakharov equations 
IH j or the equations for wave packets — the non- 
linear Schroedinger equation (NLS) and its extensions 
[3lL l32l I33L l34j . The numerical methods based on exact 
equations are quite "expensive" and thus provide a rela- 
tively low spatial resolution (typically 128 x 64, as in the 
recent work |2^| for essentially 3D waves). On the other 
hand, the applicability of the approximate equations is 
limited by the condition that the waves must not be too 
steep. To fill this gap, some new approximate, relatively 
compact explicit equations of motion for highly nonlinear 
3D waves were needed as the basis for a new numerical 
method. Recently, as an extension of the exact conformal 
2D theory a weakly 3D conformal theory has been sug- 
gested |33, which is valid for steep long-crested waves. 
Equations of this theory contain 3D corrections of the or- 
der of e = (Ixjlqj 1 ', where l x is a typical wave length, and 
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l q is a large transversal scale along the wave crests. In 
Ref. [35J, the general case was considered, with a static 
nonuniform quasi-one-dimensional bottom profile. Since 
the corresponding general equations are rather involved, 
it is natural to focus on some simple particular cases in 
the first place, for instance on the deep water limit. The 
purpose of the present work is a further development of 
the weakly 3D conformal theory and its numerical imple- 
mentation for the most simple deep water case. Our main 
result here is that we have developed a numerical scheme 
based on FFT which is sufficiently simple, accurate and 
fast simultaneously. As an application, we have simu- 
lated a freak wave formation with the final resolution up 
to 16384 x 256. 

This article is organized as follows. In Sec. II we adopt 
for the deep water limit the general weaky 3D fully non- 
linear theory described in Ref. [3!| . We also suggest some 
modification of the obtained dynamical system which re- 
sults in the correct linear dispersion relation not only 
within a small region near the fci-axis, but in the whole 
Fourier-plane. In Sec. Ill we describe our numerical 
method and present numerical results for the problem 
of a rogue wave formation. Finally, Sec. IV contains our 
summary and discussion. 

II. WEAKLY 3D NONLINEAR THEORY 

A. General remarks 

It is a well known fact that a principal difficulty in 
the 3D theory of potential water waves is the general 
impossibility to solve exactly the Laplace equation for 
the velocity potential tp(x, y, q, t), 

(fxx + Pyy + Vqq = 0, (1) 

in the flow region — 00 < y < r](x,q,t), with the given 
boundary conditions 

<P\v=t,(x,q,t) = Q, t), <p\y=-oo = 0. (2) 

(Here x and q are the horizontal Cartesian coordinates, 
y is the vertical coordinate, while the symbol z will be 
used for the complex combination z = x + iy). There- 
fore a compact expression is absent for the Hamiltonian 
functional of the system, 

H{q,i})} = ^Jdxdq J (ip 2 x + (fl + ip 2 q )dy 

— OO 

+| J r) 2 dx dq = K{rj, ^} + V{ V }, (3) 

(the sum of the kinetic energy of the fluid and the po- 
tential energy in the vertical gravitational field g). The 
Hamiltonian determines canonical equations of motion 
(see 1^3, EE EI3 , and references therein) 

m = (8H/6i/>), -in = {SH/Sri) (4) 



in accordance with the variational principle S J Cdt = 0, 
where the Lagrangian is L = J iprjt dx dq — H. 

In the traditional approach, the problem is partly 
solved by an asymptotic expansion of the kinetic energy 
K, on a sm all p arameter — the steepness of the surface 
(see Refs. [3(I|3(||38|, and references therein). As the 
result, a weakly nonlinear theory is generated, which is 
not good to describe large-amplitude steep waves [see 
Ref. [33| for a discussion about the limits of such theory; 
practically, the wave steepness ka (k is a wave number, 
a is an amplitude) should be not more than 0.1 (that 
is h/X < 0.1/ tv « 0.03)]. Below in this section we con- 
sider a different, recently developed theory j3f| (adopted 
here for the case of infinite depth) , which is based on an- 
other small parameter — the ratio of a typical length 
of the waves propagating along the ir-axis, to a large 
scale along the transversal horizontal direction, denoted 
by q [alternatively, it is the ratio of typical wave num- 
bers kq/k x in the Fourier plane (k x , k q )\. Thus, we define 
e = (lx/lq) 2 <C 1 and note: the less this parameter, the 
less our flow differs from a purely 2D flow, and the more 
accurate the equations are. A profile y = r)(x, q,t) of the 
free surface and a boundary value of the velocity poten- 
tial ip{x, q, t) = if(x, r](x, q, t), q, t) are allowed to depend 
strongly on the coordinate x, while the derivatives over 
the coordinate q will be small: \rj q \ ~ e 1 / 2 , \ip q \ ~ e 1 / 2 . 

B. Conformal variables in 3D 

In the same manner as in the exact 2D theory [Tol 
[Lj, instead of the Cartesian coordinates x and y, we use 
curvilinear conformal coordinates u and v, which make 
the free surface effectively flat: 

x + iy = z — z(u + iv, q, t), -co < v < 0, (5) 

where z(w,q,t) is an analytical on the complex variable 
w = u + iv function without any singularities in the lower 
half- plane — oo < v < 0. The profile of the free surface is 
now given in a parametric form by the formula 

X(u, q, t) + iY(u, q, t) = Z(u, q, t) = z(u + Oi, q, t). (6) 

Here the real functions X(u,q,t) and Y(u,q,t) are re- 
lated to each other by the Hilbert transform H : 

X(u,q,t) = u- HY(u,q,t). (7) 

The Hilbert operator H is diagonal in the Fourier repre- 
sentation: it multiplies the Fourier-harmonics Yk{q,t) = 
J Y(u, q, t)e~ lku du by i sign k, so that 

HY(u,q,t) = J[isignk}Y k (q,t)e lku (dk/27v). (8) 

The boundary value of the velocity potential is tp\ v =o = 
i/j(u,q,t). 

For equations to be shorter, below we do not indicate 
the arguments (u, q, t) of the functions if), Z and Z (the 
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overline denotes complex conjugate). The Lagrangian for 
3D deep water waves in terms of the variables ip, Z, and 
Z can be written as follows (compare with |3fij): 



ZtZ u — ZtZ u 



2 1 



21 



2 r 



ip dudq — K,{ip, Z, Z} 
a. + Z.. 



du dq 



A 



H 



Z-Z 
2i 



Z + Z 



dudq, (9) 



where the indefinite real Lagrangian multiplier A(u,q,t) 
has been introduced in order to take into account the 
relation JJJ . Equations of motion follow from the varia- 
tional principle 5 A = 0, with the action A = J Cdt. So, 
variation by Sip gives us the first equation of motion — 
the kinematic condition on the free surface 



lm{Z t Z u ) = {SK/Sip). 



(10) 



Let us divide this equation by \Z U \ 2 and use analytical 
properties of the function Z t /Z u . As a result, we obtain 
the time-derivative-resolved equation 



Z t = iZ u (l + iH) 



(SIC/Sip) 



(11) 



Further, variation of the action A by SZ gives us the 
second equation of motion: 



Tp u Z t - 1p t Z u 



2 1 



(12) 



After multiplying Ea. (jT2*|) by — 2iZ u we have 
[ip t + glmZ]\Z u \ 2 - ip u Z t Z u = 



i-#)A-2t (jP)Z U) 
. . (13) 

where A is another real function. Taking the imaginary 
part of Ea. l(T3"|) and using Ea. lfTTTjl . we find A: 



A 



" , SIC' 








+ 2 Re 


[(S)-] 



(14) 



After that, the real part of Ea. (|13|) gives us the Bernoulli 
equation in a general form: 



tpt = -g Im Z - ip u H 
Im ( (1 



(6tC /Sip) 



Z u \ 2 

H) [2(SKjSZ)Z u + (5tC/8ip)ip u f) 

\zj 2 



(.15) 



Equations (|llfl and (|15|l completely determine evolu- 
tion of the system, provided the kinetic energy functional 
lC{ip, Z, Z} is explicitly given. It should be emphasized 
that in our description a general expression for K, re- 
mains unknown. However, under the conditions \z q \ <C 1, 



\ip q \ <C 1, the potential ip(u, v, q, t) is efficiently expanded 
into a series on the powers of the small parameter e: 



p=p^ +pW + 



(pW ~ e r ' 



(16) 



where p( n+lS> can be calculated from (p^ , and the zeroth- 
order term ipW — Re <p(w, q, t) is the real part of an easily 
represented (in integral form) analytical function with 
the boundary condition Re^|„ = o = ip(v>,q,t). Corre- 
spondingly, the kinetic energy functional will be written 
in the form 

K = /C (0) + /C (1) + . . . , /C (,l) ~e™, (17) 
where K.^°'{ip} is the kinetic energy of a purely 2D flow, 

IC (0) m = \f[{p^) 2 + {p { v 0) f]dudvdq 



ipHip u du dq, 



(18) 



and the other terms are corrections due to gradients along 
q. Now we are going to calculate a first-order correction 



C. First-order corrections 

As a result of the conformal change of two variables, 
the kinetic energy functional is determined by the expres- 
sion 

£ = \ J H + + AQ ■ V^) 2 ] du dv dq, (19) 

where the Cauchy-Riemann conditions x u = y v , x v = 
—y u have been taken into account, and the following no- 
tations are used: 

J ee \z u \ 2 , (Q • Vip) = atp u + bcp v + <p q , 



Consequently, the Laplace equation in the new coordi- 
nates takes the form 

Puu i Pvv 

V ■ (QJ(Q • Vp)) = 0, (20) 

with the boundary conditions 

tp\ v=0 = ip{u,q,t), tp\ v =-oa=0. (21) 

In the limit e <C 1 it is possible to write the solution as 
the series <|16[) . with the zeroth-order term satisfying the 

2D Laplace equation p^u + Pvv = 0, p\fl = ip(u,q,t). 
Thus, it can be represented as p^ = Re (p{w, q, t), where 

r+°° fju 
<P{w,q,t) = (1- sign k)M^t)e tkw — , (22) 

J-oo 27r 
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and ipk(q,t) = J ip(u,q,t)e %ku du. On the free surface 

4>{u + iO, q, t) ee q, t) = (l + iH)ip(u, q, t). (23) 
For all the other terms in Eci. (|16J) we have the relations 
V-(QJ(Q-Vv> (n) )) = (24) 



Jn+l) + (n+l) 
tuu 1 tvv 



and the boundary conditions tp 



(n+l) 



v=0 



0. Noting 



that Jipu^Pu^ + ifv <Pv) dudv = (it is easily seen 
without explicit calculation of (p^' after integration by 
parts), we have in the first approximation 

/C« = \ / 'j(^°) + mpW + bp^) 2 dudvdq 



Re ^ 




<i>uZ q \ 









du dv dq. (25) 



Since z(w) and <j)(w) are represented as z(it + iv) = 
e~ kv Z(u) and <f>(u + iv) = e v ^(u), we can use for v 
integration the following auxiliary formulas: 

pO rO 

du 



[e- kv A{u))[e- kv B{u)] dv= - 



A k B k dk 



B(u)d- 1 A(u))du 



2k 2tt 
(26) 



Now we apply the above formulas to appropriately de- 
composed Eq.J35j) and, as a result, we obtain an expres- 
sion of a form K.^ — Z, Z}, where the functional 
T is defined as follows (compare with |35jl 



T = g / (^«* g - Z q ^ u )d-\Z u ^ q - Z q V u )dudq 



y/zj (z-u) 



-{Z - u) [{Z u V q - Z q V u ) 2 /Z u ]j dudq.{27) 

Here the combination (Z — u) is written instead of Z 
just for convenience, as it is finite at the infinity. Ac- 
tually the equations of motion "do not feel" this dif- 
ference. In the derivation we have used the identity 
J[(Z u ^ q — Z q ^ u ) 2 / Z u ]udu = 0, which holds because the 
integrand is analytical in the lower half-plane. From here 
one can express the variational derivatives (SK.^ /Sip) 
and (SIC^ /SZ) by the formulas 



dtp 



= 2 Re 



(1-iH) 



5JF_ 
S^ 



SJC^_ 

SZ 



ST_ 
SZ' 



(28) 



The derivatives (SJ-/8^) and (5J-/8Z) are calculated in 
a standard manner: 

S* ~ 8 q 



{z u m q - z q m u ) 



+d u [{V q - Z q V u /Z u ){Z-u)] 

l - z u d q la-^z^g-Zg^) 



+ (* q -Z q * u /Z u )(Z-u) , (29) 



ST_ 

Jz 



(Z u V q - Z q V u ) 



+d u [(* q - z q * u /z u )(z-u)} 



d~'{z u m q - z q ^ u ) 



i 

To 



+(V q -Z q $ u /Z u ){Z-u) 
d u [(V q - Z q y u /Z u ) 2 (Z-u)} 



-(* q -Z g * u /Z u )*Z v 



(30) 



Now one can substitute (SIC/ Sip) « -Htp u + (SK^/Sip) 
and (SK./SZ) m (SK.^ /SZ) into the equations of motion 
(| 1 1 ft and H15JI . Thus, the required weakly 3D equations 
for deep water waves are completely derived: 

Z = u+(i-H)Y(u,q,t), Z u = 1 + (i - H)Y U . (31) 
Z t = -iZ u (l + %A) \[Hip u - {SJC^/S^)]/\Z U \ 2 ] , (32) 



ipt = -gY+ ip u H 



[Hip u -{5K,V/5i>)]/\Z u \ 



H 



il> u [H^-(6K.W/6il>)] /\Z U 



2 Re ((H + i)[Z u {5JC^/5Z)])/\Z u 



(33) 



D. Modification of the Hamiltonian 



It can be easily obtained that the linear dispersion re- 
lation for the system (|27(l - (|33|l is 



u)\k,m) = g\k\ 1 + 



1 

2~P 



(34) 



where m is the wave-number in the g-direction (the wave- 
number k in the u-direction was introduced earlier) . Ob- 
viously, here we have the first two terms from the expan- 
sion of the exact 3D linear dispersion relation 



u) 2 (k, m) — g\/k 2 + m 2 



(35) 



on the powers oim 2 /k 2 < 1. Thus, the system (j2Tf - (pf5)) 
has a non-physical singularity in the dispersion relation 
near the m-axis in the Fourier plane (fc, m), where k/m — > 
0. Therefore, for convenience of the numerical modeling, 
some regularization could be useful in the approximate 
Hamiltonian K.^ + + V, where V is the potential 
energy, 



p = £ 

2 



z-z 



21 



Z u + Z u 



du dq. 



(36) 



5 



Also, the correct linear dispersion relation is strongly de- 
sirable, since it determines which waves interact by non- 
linear resonances. Of course, a possible regularization is 
not unique as we keep only zeroth- and first-order terms 
on e = (m/fc) 2 in the Hamiltonian. Below we suggest a 
modification which adds terms of the order C(e 2 ) to the 
approximate Hamiltonian V + JC^ + K,^ (a modified 
Hamiltonian 7i = V + K, will remain valid up to e) . First 
of all, instead of the the functional K.^ we use another 
functional, K.^ — > (1/2) J tpG^xp dudq, where the linear 
operator Go is diagonal in the Fourier representation: 



G (k,m) = VP + m 2 -!-^ 



2 fc 2 + TO 2 ' 



(37) 



Besides that, we change the operator id u 1 = 1/fc in the 



first line of Eq. (|T7Jl by the operator fc/(fc 2 



iduA^ 1 , which is less singular. As the result, we have 
the modified approximate kinetic-energy functional in 
the form 



*"5 



tpGoip du dq + 



(38) 



T = - \ (Z u V q - Z q ^ u )d u A^{Z u ^ q - Z q y u )dudq 



+ j [(z u v q -z q y u ) 2 /z u ] (z- 



-{Z-u) [(Z u * q - Z q ^ u ) 2 /Z u ]jdudq. (39) 

The linear dispersion relation resulting from K, is correct 
in the whole Fourier plane. Besides that, the zeroth- 
and the first-order terms on e in K, are the same as in 
K.^ + KS 1 '. The system of equations, corresponding to 
the modified Hamiltonian H, consists of Eq. (|31J) and the 
following equations: 



Z t = iZ u (l + iH) (8K/6il;)/\Z 



, I 2 

J U\ ! 



(40) 



Jz = 



d u Az\Z u y q - Z q V u ) 



-^> q -Z q ^> u /Z u ){Z-u) 



v u d q d u A^(z u ^> q - z q ^ u ) 



-(* q -Z q *jZ u )(Z-u) 



I 

16 



d u [{* q - z q * u /z u y{z - u)] 



-(V q -Z q V u /Z u ) 2 Z u . (44) 



These equations were integrated numerically as de- 
scribed in the following section. It is interesting to note 
that we also considered another choice for the regulariza- 



= tion, with G = (fc 



V2)(fc 2 



2 )-V2 i ns tead of \k\ 



in K.(°>, and (k 2 + m 2 ) -1 / 2 instead of l/\k\ in the first 
term in K,^ . The numerical results were found very close 
in both cases, since the wave spectra were concentrated 
in the region m 2 /k 2 <C 1. 



III. NUMERICAL METHOD AND RESULTS 

In our numerical simulations, we used the following 
procedure. A rectangular domain L u x L q with the peri- 
odic boundary conditions in the (u, g)-plane was reduced 
to the standard dimensionless size 2ir x 2ir [the aspect 
ratio £ = (L u /L q ) 2 was taken into account]. This stan- 
dard square was discretized by N x L points (u n ,qi) — 
2n(n/N,l/L), with integer n and /. The time variable 
was rescaled to give g — 1. As the primary dynami- 
cal variables, the (complex) Fourier components Yfc m and 
Pkm = (k 2 +sm 2 ) 1 ^ 4 ipkm were used, where k and m were 
integer numbers in the limits < k < K, —M < m < M , 
with K w 3A/8, M « 3i/8. For ne gativ e fc, the prop- 
erties F_ fe) _ m = Y km and P-k-m = Pkm were implied. 
It should be noted that after rescaling u, q, and t, the 
linear dispersion relation has been Ukm — (fc 2 + em 2 ) 1 / 4 , 
and that the combinations 



ipt = lm(p.-iit)[2(6F/8Z)Z» + (8K/54)il> u ])/\Z u \ 



5JF_ 



gY-^ u H {6K/6$)/\Z. 



2 

u I ? 



^ = G </> + 2Re 
dip 
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d u A^{Z u ^ q -Z q ^ u ) 



+(V q -Z q V u /Z u )(Z-u) 



' 77 Z u d„ 



d u A^(Z u ^ q - Z q ^ u ) 



+ (V q -Z q * u /Z u )(Z-u) , (43) 



(41) 



(42) 



km 



iP, 



km 



\j2uJkr, 



in the linear limit coincide with the normal complex vari- 
ables [HIH. 

Using the FFTW library 0, the quantities Z(u n , q t ), 
^(u n ,qi), their corresponding u- and ^-derivatives, and 
Goip(u n , qi) are represented by 2-dimcnsional complex ar- 
rays. The variational derivatives (|43|l . Q44JI. 142|l and 
the right-hand-sides of Eqs. (|40|l - (|41|) are then computed 
and transformed into Fourier space in order to evalu- 
ate the multi-dimensional function which determines the 
time derivatives of Y km and Pkm ■ This function is used in 
the standard fourth-order Runge-Kutta procedure (RK4) 
for the time integration of the system. After each RK4 
step, a large- wave-number filtering of the arrays Yk m and 
Pkm is carried out, so that only the Fourier components 
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FIG. 1: (Color online) [A]. Map of the free surface at t — 0. 



with < k < K e ff, —M e ff < m < M e ff are kept, where 
K e ff w N/4 and M e // ~ L/4. To estimate the accuracy 
of the computations, conservation of the total energy Ti. 
and the mean surface elevation are monitored. 

Typically at t = we put N = 2 12 , L = 2 6 , and 
the initial time step t = 0.01. During the computa- 
tion, as wave crests become more sharp and the spec- 
tra get broadened, we several times adaptively double 
N (together with K, K e ff) and L (together with M, 
M e ff), with the time step half-decreasing when N is dou- 
bled. At the end, when a giant wave is formed, we have 
N = 2 13-14 , L = 2 7 ' 8 . As a result of such adaptive 
scheme, the conservation of the total energy is kept up 
to 5-6 decimal digits during most part of the evolution. 
Only at a very late stage, when N and L are not allowed 
to double anymore, the conservation is just up to 3-4 
digits, and the filtering of the higher harmonics becomes 
more influential. 

The efficiency of the above described numerical method 
crucially depends on the speed of the FFT routine, since 
most of the computational time (approximately 80%) is 
spent in the Fourier transforms. To compute the evolu- 
tion through the unit (dimensionless) time period, on a 
modern PC (Intel Pentium IV 3.2 GHz) it takes three- 
four minutes with N = 4096, L = 64, r = 0.01, and 
more than one hour with N = 8192, L = 256, r = 0.005. 
The total time needed for a single numerical rogue-wave 
experiment is about 3-5 days. 

Two numerical experiments are reported below for hor- 
izontal physical box sizes of 4 x 4 km (the final resolution 
was 16384 x 256 in both). In the first of them, the main 
mechanism for a big wave formation is the linear disper- 
sion, when a group of longer waves overtakes a group of 
shorter waves, and their amplitudes are added (see Q 
for a discussion). However, due to nonlinear effects, in 
maximum the crest is higher than simply the sum of two 
amplitudes. In the second experiment, a giant wave is 
formed by an essentially nonlinear mechanism (due to 
Benjamin-Feir instability) from a slightly modulated pe- 
riodic wave, close to a stationary Stokes wave. 
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FIG. 2: [A]. Initial wave profiles for eight equidistant values 
? = L 9 (0,l/8,...,7/8). 
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FIG. 3: (Color online) [A]. Map of the free surface at t = 59.2 
(the physical time is 7 min 57 sec). 
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FIG. 4: [A]. Wave profiles at t = 59.2. 
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FIG. 7: (Color online) [B]. 
172.8. 



Map of the free surface at t 



FIG. 6: [B]. Rogue wave at t = 172.8 (23 min 13 sec) 



A. "Linear" big wave 

In the first numerical experiment (referred to as [A]), 
the initial state was a composition of two spatially sepa- 
rated wave groups, as shown in Fig^ Two typical dimen- 
sionless wave numbers are present, ki — 17, correspond- 
ing to the wave length Ai « 235 m (at the central part of 
Fig.UJ and fc 2 = 23 (A 2 « 174 m). At t = Q the highest 
crests were about 5 m in both groups (see Fig. 01. Due 
to difference in the group velocities (c gr ~ |fc| -1 / 2 ), the 
longer waves move faster and after some time overtake 
the shorter waves. As a result of almost linear superpo- 
sition, big waves with high crests and deep troughs are 
formed, separated by the spatial period 27r/|fc 2 — with 
the maximum amplitude about 12 m, which is approxi- 
mately equal to the sum of the individual amplitudes (the 
nonlincarity slightly increases the maximum height). The 
corresponding numerical results are presented in Figs|21 
andgl 



B. Nonlinear big wave 

Our second numerical experiment (referred to as [B]) 
is a weakly 3D analog of the 2D numerical experiments 
performed by Zakharov and co-workers 0, U, when a 
slightly modulated periodic Stokes wave evolves to pro- 



duce a giant wave. However, with two horizontal dimen- 
sions we could not achieve the same very high resolution 
as Zakharov with co-workers did for the case of a single 
horizontal dimension (16384 x 256 in our experiment [B] 
versus 2 • 10 6 x 1 in Ref. Q). Our computations were 
terminated well before the moment when a giant wave 
reached its maximum height and began break, since the 
number of the employed Fourier modes became inade- 
quate. The accuracy could be better with larger N and L, 
and with a smaller r, but it required much more memory 
and computational time. In general, to resolve in confor- 
mal variables a sharp wave crest with a minimal radius 
of curvature p and with the asymptotic angle 27r/3 (as 
in the limiting Stokes wave), the required number K e ff 
should adaptively vary as p _3 / 2 A 1 / 2 . The power —3/2 
results in strong difficulties when p is small. Another im- 
portant point is that Zakharov and co-workers were able 
to re-formulate the purely 2D equations in terms of the 
"optimal" complex variables R = 1/Z U and U = i^ u /Z U: 
thus obtaining very elegant and compact cubic evolu- 
tion equations (the Dyachenko equations [13 ). In our 
3D case, a similar simplification seems to be impossible, 
and we dealt directly with the original conformal vari- 
ables Z and ^ . Nevertheless, our results are sufficiently 
accurate to reproduce the fact of a giant wave formation. 

As the initial state for the experiment [B], we took 
a weakly modulated periodic wave with the main wave 
number k = 19 (A ~ 210 m) and with a few first har- 
monics on 2k, 3fc, etc., similar to a Stokes wave (not 
shown). After some period of evolution, the Benjamin- 
Feir instability developed and resulted in formation of a 
big wave, with the amplitude 13 m at t — 172.8 versus 
the initial maximum amplitude 5 m (see Figs. 15171 and 
compare with Ref. (3J)- The peak-to-trough height h* 
of this computed rogue wave was approximately 20 m 
at t = 172.8 (the steepness parameter h*/\ sa 0.1), and 
it was still growing at that moment (so, at t = 173.0 
we observed the amplitude 14 m, but the accuracy was 
already not sufficient). It is interesting that this numeri- 
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cal solution has a well denned envelope until a very final 
stage of evolution. Thus, our equations may serve to test 
the simplified wave-packet models like the extended NLS 
equations pul 132. 1331 ] . 

IV. SUMMARY AND DISCUSSION 

We have developed an efficient numerical method for 
modeling the rogue wave phenomenon. The underlying 
theory for the method is the weakly 3D formulation of 
the free surface dynamics reported earlier |35j ] with slight 
modifications. In particular, the Hamiltonian has been 
regularized in a way to give the exact linear 2D disper- 
sion relation in the entire Fourier space, and a filter re- 
moving shortest waves has been added in the numerical 
implementation. With these techniques, weakly three- 



dimensional effects could be included in simulations of 
rogue wave formation as illustrated in the two examples 
given. In particular, the genuinely non-linear 2D insta- 
bility reported by Zakharov et al. 0, could be verified 
in the weakly 3D regime despite the relatively low reso- 
lution (compared to 2D) of 16384 x 256 points used here. 

The results indicate that the assumption of weak vari- 
ation in the third direction holds even in the late stage 
of rogue wave formation, which demonstrates the consis- 
tency of the expansion in e = (l x /lq) 2 and thereby the 
applicability of the present theory. 

Planned further steps in the continuation of this 
work are a more efficient computational implementation 
through parallelization of the code and the inclusion of 
additional effects like a bottom profile, which is already 
covered by the formalism reported earlier |35| . 
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